High Accuracy Many-Body Calculational Approaches for Excitations in Molecules 
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Two state-of-the-art computational approaches: quantum Monte Carlo (QMC), based on accurate 
total energies, and GW with exciton effects (GW-BSE) , based on perturbation theory are employed 
to calculate ionization potentials, electron affinities, and first excited singlet and triplet energies for 
the siiane and methane mofecufes. Resufts are in excellent agreement between these dramatically 
different approaches and with available experiment. The optically forbidden triplet excitation in 
siiane is predicted to lie roughly f eV higher than previously reported. For methane, the impact of 
geometry relaxation is shown to be 2 eV for excited states. Further, in the GW-BSE method, we 
demonstrate that inclusion of off-diagonal matrix elements in the self-energy operator is crucial for 
an accurate picture. 



Optical excitations play a fundamental role in techno- 
logical applications such as dye chemicals, photovoltaics, 
laser technology, and catalysts for chemical reactions. 
An accurate evaluation of excited states requires more 
sophisticated approaches than simply a straightforward 
application of the mainstream ground state methods such 
as density functional theory and Hartree-Fock. In addi- 
tion, many traditional approaches for calculating excited 
states can be limited by a range of factors. For exam- 
ple, configuration interaction with single excitations is 
missing correlation effects and is often inaccurate, optical 
transitions modeled as quasiparticle energies corrected by 
excitonic effects rely on an approximate description of the 
electron-hole interaction, and basis set correlated wave 
function approaches are limited by their scaling with the 
number of correlated electrons. 

Two promising new approaches have been developed 
to calculate accurate excitation energies: (1) quantum 
Monte Carlo (QMG), based on a stochastic solution 
to the many-electron Schrodinger equation, can provide 
chemical accuracy (i.e., 0.05 eV) for the total energy 
while at the same time scales as where N is the number 
of valence electrons (see, e.g., [0,^), and (2) many-body 
Green's function perturbation techniques that employ 
the GW approximation |^,^ for the electron self-energy 
operator to calculate the single-particle excitation spec- 
trum, followed by a solution to a Bethe-Salpeter equation 
(BSE) for electron- hole excitations 

In this Letter, we demonstrate excellent agreement be- 
tween QMC, GW-BSE, and available experiment for ex- 
citation energies of siiane (Sill4) and methane (CII4). 
Calculated ionization potentials, electron affinities, and 
first excited singlet and triplet energies also agree very 
well with highly reliable standard quantum chemistry 
approaches; however, the major advantage of our two 
current approaches is numerical efficiency. Both meth- 
ods require moderate computation times, and their scal- 
ing with system size is to our knowledge superior to all 



other quantum-chemistry approaches for computing ex- 
cited states. The excellent agreement between QMC and 
GW-BSE is encouraging given the dramatically different 
methodologies involved (i.e., the former relies on accu- 
rate total energy differences while the latter is based on 
perturbation theory). 

In our QMG approach (7 - 11 1 , variational Monte Carlo 
is employed to find an optimized correlated many-body 
trial function. This trial function, ^t{R), is a product of 
Slater determinants, and a correlation factor p2]. 



= X! d^^li^n exp 
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where / corresponds to the ions, i,j to the electrons 
and Til, Tji, Tij to the relevant separations. Parametriza- 
tion and optimization of u(rii, rji, rij), which represents 
the electron-electron and electron-electron-ion correla- 
tions, is described in Ref . jl^ . In the Slater determinant 
part, we employ natural orbitals rather than Hartree- 
Fock (HF) or density functional orbitals . To eliminate 
most of the remaining variational bias we use the fixed 
node diffusion Monte Carlo method, which is based on 
the property that the operator exp(— ri?), where H is the 
Hamiltonian, projects out the ground state of any trial 
function with the same symmetry and non-zero overlap, 
subject to the constraint that the nodes are unchanged 

All QMC results presented here are from the diffu- 
sion Monte Carlo approach. 

The GW-BSE Green's function perturbation approach 
consists of three successive steps. In the first step a lo- 
cal density approximation (LDA) | [T5[ | calculation is per- 
formed to obtain the electronic ground state. In the sec- 
ond step the LDA results are used to construct the elec- 
tronic self-energy operator S (also known as the mass 
operator). This is done within the GW approximation 

where G is the (LDA) single-particle Green's func- 
tion and W is the screened Coulomb interaction based on 
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RPA dielectric screening. Solving the equation-of-motion 
for the single-particle Green's function yields the quasi- 
particle states (i.e., the occupied (hole) states v and the 
empty (electron) states c) of the molecule. In the third 
step the electron-hole interaction {vc\K'^'^\v' c') between 
the occupied states v and the empty states c is evaluated. 
The coupled electron-hole excitations |S') = X]«c ^^d'^^c) 
and their excitation energies Sl"^ then result from the 
Bethe-Salpeter equation 
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For details, we refer to Refs. [|6|,|l6[. 

We have also calculated excited state energies as total 
energy differences using LDA and the generalized gra- 
dient approximation (GGA) with the PW91 functional 
[0. In addition, the small number of correlated elec- 
trons in the present systems permits us to employ sev- 
eral of the most accurate quantum chemical approaches: 
coupled cluster with single, double, and perturbationally 
triple (CCSD(T)) excitations ||l^, and complete active 
space self consistent field (CASSCF) with an active space 
of 8 electrons and 9 active orbitals. For these methods, 
extensive basis set tests showed that energy differences 
are converged to better than 0.02 eV with the cc-pVQZ 
p9[ basis set. All calculations presented here other than 
QMC and GW-BSE have been carried out using the G98 
package and the cc-pVQZ basis set. 

It is important to distinguish between the two different 
methodologies employed here. Using total energy meth- 
ods, the ionization potential (IP) and electron affinity 
(EA) are computed as total energy differences. In the 
Green's function, the QP energies are obtained in a dif- 
ferent manner, namely as a solution of a quasiparticle 
Dyson's equation 

The single-particle HF and LDA energy spectra for 
SiH4 are shown in Fig. |^. While both methods are in 
qualitative agreement with one another, there is a large 
quantitative discrepancy, as expected: HF usually over- 
estimates HOMO-LUMO gaps and eigenvalue spacings, 
whereas LDA typically underestimates these quantities. 
Note that since the first excitations are from t2 to ai or- 
bitals, the ground and first excited singlet states of SiH4 
are of different symmetry. 

It is reassuring that all of the methods (see Table |) 
are in good agreement with experiment for the ionization 
potential of SiH4. The LDA, GGA, and HF results are 
~ 0.5 eV lower than experiment, indicating a slight bias 
towards the ionized state; however, note the dramatic im- 
provement when compared to their single-particle highest 
occupied eigenvalue. The calculations indicate a small, 
negative (i.e., nonphysical) electron affinity ranging from 
-0.1 to -0.6 eV. 

Within the GW-BSE approach, the behavior of the ai 



LUMO state deserves careful discussion (see Table ||). 
In the underlying LDA calculation, the LUMO energy is 
below the vacuum level (see Fig. 1) and is thus given 
by a localized, bound wave function. Within first-order 
perturbation treatment of the self-energy operator E, the 
orbitals are kept as LDA orbitals and the LUMO energy 
is shifted above the vacuum level by 1.7 eV (see second 
column of Tab. ||). Being unbound now, the QP state 
changes its wave function from a localized to a more ex- 
tended orbital. This is described by going beyond the 
diagonal evaluation of S in the LDA basis in the GW 
calculation and taking off-diagonal elements of S into 
account, which are in fact non-zero. This effect, which 
was not accounted for in a previous work lowers the 
LUMO state by as much as 0.8 eV (see third column of 
Tab. p. The LUMO energy of 0.3 eV is now in good 
agreement with the results of the other methods. Note 
that the LUMO state (being unbound) converges slowly 
with respect to the basis, due to its delocalized nature. 
The low-energy electron-hole excited states, on the other 
hand, that are the focus of this paper, are spatially lo- 
calized since the electron is attracted to the hole; these 
excited states are much easier to converge. 

Our results for SiH4 excitation energies for the lowest- 
energy triplet and singlet excited states are given in Ta- 
ble III. Calculation of spin-singlet excitations by total 
energy differences poses a serious challenge to traditional 
approaches. For example, the first excited singlet state 
is described by a 2 determinant wavefunction. Accurate 
correlation of this state was therefore only possible using 
the CASSCF method. The triplet state is significantly 
easier to converge within traditional approaches owing 
to the fact that it is the lowest energy state of the given 
spin multiplicity. The QMC approach does not encounter 
the same difficulties as the traditional approaches for the 
singlet state since the spin configuration is given by spec- 
ification of the slater determinant ~ singlet and triplet 
states are equally as simple to simulate, provided each 
wavefunction is orthogonal to the ground state. In the 
case of the excited singlet state, the DMC energy was 
found to be rather insensitive (i.e., <0.1 eV) to the num- 
ber of determinants in the trial function, indicating that 
the nodes of a single determinant closely resemble those 
of the 2-determinant wavefunction which has the correct 
spatial symmetry. 

The SiH4 spin singlet results are in good agreement 
with experiment |21 for the three theoretical methods 
listed in Table 



III 



The optically forbidden triplet ex- 
citation is difhcult to obtain from experiments. Using 
electron collision spectroscopy, slight spectral structures 
were found at around 7-8 eV We believe that the 

agreement among the theoretical methods presented here 
strongly suggests that the triplet excitation lies roughly 
1 eV higher than previously reported. 

The GW-BSE results of Table |!| differ from the ones 
of Ref. [| ( see also Table ||). The reason is again due 
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to the off-diagonal matrix elements of the self-energy op- 
erator when expanded in the LDA orbitals, which had 
been neglected in Ref. As discussed above, the full 
self-energy operator allows the unoccupied states to be- 
come more delocalized. Therefore the overlap between 
the hole and electron states, as well as the electron-hole 
interaction, is significantly reduced. Concomitantly, the 
excitation energies (column 3 in Tab. ||) are not re- 
duced compared to column 2 (as the lower LUMO energy 
suggests), but are in fact increased. It is interesting to 
note that the triplet excitation energy is more affected 
than the singlet, which reduces the singlet-triplet split- 
ting from 1.1 eV to 0.7 eV within GW-BSE. 

Methane and silane are isoelectronic and similar in 
their valence shell structure; however, the stronger po- 
tential of carbon compared to silicon causes wavefunc- 
tions and densities to be less smooth which in turn can 
lead to difficulties and inaccuracies within ab initio cal- 
culations. We have carried out the same calculations de- 
scribed above for CH4. Our results, shown in Table IV, 



again demonstrate strong agreement between QMC and 
GW-BSE. 

An important issue in the discussion of molecular spec- 
tra is the interaction between electronic excitations and 
the molecular geometry. Despite the high dimensionality 
of a five-atom molecule, we are able to discuss the most 
important aspects of geometric relaxation in the excited 
states and their consequences for the total energy. Table 
IV illustrates the impact of relaxation from the ground 
state Td symmetry for each excited state molecule. Note 
that the optimization has been restricted to sym- 
metry. Test calculations show that further lowering the 
symmetry results in a total energy reduction no larger 
than ~ 0.1 eV. 

For all three excited states shown here, geometry re- 
laxation yields an enormous lowering of the total energy 
by as much as 1.9 eV. A minor contribution (about 0.2- 
0.4 eV) is related to the increase of the C-H distance in 
the excited states. The main contribution results from 
breaking the symmetry of the ground state. This ef- 
fect can be understood from the molecular orbital scheme 
(see Fig. |^). Within the ground state, all three degener- 
ate HOMO levels are fully occupied, i.e. the ground state 
is non-degenerate. In the excited states, one of the three 
HOMO states becomes half-filled, i.e. the excited states 
are three-fold degenerate. The symmetry-breaking from 
Td to C2V lifts the degeneracy of the three HOMO levels. 
In the single-particle spectrum of Fig. |l] two HOMO lev- 
els are shifted down and one is shifted up, forming the 
new HOMO level of the molecule in the reduced sym- 
metry. In the optimum geometry of the excited states 
this level-splitting is about 4.2 eV, although the center- 
of-mass of the three levels remains roughly the same. The 
new HOMO level, which is half-filled in all three excited 
states, is non-degenerate, so the excited states are non- 
degenerate as well. 



A closer analysis of the excited-state geometry shows 
that the increase of the C-H bond length is relatively 
small (no more than 0.07 A). The main effect results 
indeed from the Jahn- Teller distortion of the molecule, 
which changes the bond angles from 109.47 degree (tetra- 
hedral angle) to 94 degree for the small and 118 degree 
for the large bond angles. This distortion increases the 
total energy of the ground state by 1.5 eV. In the ex- 
cited states, however, this increase of the total energy 
is overcompensated by the splitting of the three HOMO 
levels described above, i.e. by the upwards shift of the 
new HOMO level and a correspondingly large reduction 
of the excitation energies. In the experimental spectrum 
for CH4 H^, a weak shoulder is observed at about 10 eV, 
corresponding to a vertical excitation. Compared with 
the onset at 8.52 eV, this indicates a relaxation gain of 
about 1.5 eV, supporting our calculated results. 

Similar effects on the excitation energy are also found 
for SiH4. The gain in total energy is 0.9 eV, i.e. slightly 
smaller than in GH4. In the measured absorption spec- 
trum this corresponds to the difference between the en- 
ergy of maximum absorption and the low-energy onset. 
This difference amounts to 0.6 eV in the absorption mea- 
surements of Ref. pl| . Our minimum-energy transition 
for SiH4 is 8.3 eV (in GW-BSE), in excellent agreement 
with the measured onset of the spectrum at 8.2 eV. 
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TABLE II. Quasiparticle levels and lowest electron-hole 
excitation energies (eV), calculated for SiH4 within LDA and 
within two evaluations of the GW-BSE approach. Unlike the 
second column, the GW-BSE results of the third column in- 
cludes off-diagonal matrix elements of the self-energy operator 
and allows for changes in the one-particle wave functions. 



LDA GW-BSE (diagonal S) GW-BSE (full S) 



HOMO energy -8.4 -12.7 -12.7 

LUMO energy -0.6 1.1 0.3 

Es 7.8 8.8 9.2 

Et 7.8 7.7 8.5 



TABLE HI. Neutral excitation energies (eV) for the first 
excited singlet and triplet states of SiH4 computed as 
E(excited state) - E(ground state) for all methods except 
GW-BSE. A corresponds to the singlet-triplet splitting. The 
experimental result denotes the energy of maximum ab- 
sorption for the state studied here, which corresponds to ver- 
tical excitation. 
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FIG. 1. Single-particle Hartree-Fock and local density ap- 
proximation eigenvalue spectra (eV) for the SiH4 molecule. 
Note that the CH4 molecule has the same orbital structure. 



TABLE I. Ionization potential (IP) and quasiparticle (QP) 
gap (eV) for SiH4, computed as E(N-1)-E(N) and E(N-f 1) + 
E(N-l) - 2E(N), respectively, for all methods except GW. 



IP QP Gap 



HF 12.3 12.6 

LDA 12.1 12.4 

BPW91 12.1 12.5 

CCSD(T) 12.7 13.3 

GW 12.7 13.0 

QMC 12.6(1) 12.8(1) 

EXP" 12.6 — 



TABLE IV. Ionization potentials (IP) and energies of the 
lowest singlet (Es) and triplet (Et) excited states for CH4 in 
both the ground state (Td) and relaxed {C2v) symmetries. All 
energies are in eV and are given with respect to the ground 
state. The experimental result [^ denotes the onset of the 
spectrum, which corresponds to the minimum energy transi- 
tion including structural relaxation. 
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